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Abstract 



Lattice Boltzmann simulations are used to investigate spinodal decompo- 
sition in a two-dimensional binary fluid with equilibrium lamellar and droplet 
phases. We emphasise the importance of hydrodynamic flow to the phase sep- 
aration kinetics. For mixtures slightly asymmetric in composition the fluid 
phase separates into bulk and lamellar phases with the lamellae forming dis- 
tinctive spiral structures to minimise their elastic energy. 
PACS numbers: 47.11. +j; 64.75.+g 
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I. INTRODUCTION 



Complex fluids such as microemulsions, foams, and colloidal suspensions provide a wealth 
of interesting physical phenomena. Their use is ubiquitous in the processing, energy and 
chemical industries. Therefore it is important to explore ways in which their properties can 
be modelled numerically to test physical theories, predict behaviour and provide feedback 
to industrial planning. 

The modelling of complex fluids is not an easy task as both the rheology and the compli- 
cated phase behaviour of the fluid must be incorporated. Molecular dynamics simulations, 
whilst providing an accurate picture of the microscopic physics, are usually too computer- 
intensive to address hydrodynamic time scales. In computational fluid dynamic solutions of 
the Navier-Stokes equations the specific behaviour of a given fluid can only be input via ad 
hoc constitutive relations. 

Recently new methods of simulating fluid flow have become available. These in- 
clude lattice-gas cellular automata^, dissipative particle dynamicsi, and lattice Boltzmann 
simulations^. The aim is to reproduce the physics of fluid flow, primarily mass and mo- 
mentum conservation, while including the important features of the underlying microscopic 
physics. These approaches have been successfully applied to several complex fluids including 
polymer solutions!, particulate suspensions!! and microemulsionsi. However application of 
the techniques to model complex fluids is still in its infancy and validation of, and comparison 
between, the different methods is still needed. 

In this paper we concentrate on lattice Boltzmann simulations. A fluid is modelled on 
a mesoscopic length scale by means of distribution functions which evolve according to a 
discretised version of a simplified Boltzmann equation. The correct equilibrium behaviour 
is imposed by inputting the correct thermodynamics such that the system evolves to the 
minimum of a chosen input free energy!'!. 

We use the lattice Boltzmann approach to explore the kinetics of phase separation of 
a model, two-dimensional, fluid with lamellar and droplet equilibrium phases. Our main 
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conclusions are: 

1. When the system is quenched to the lamellar phase a glassy metastable state of tangled 
lamellae forms. However, hydrodynamic flow alleviates the frustration and allows the 
system to attain the striped ground state. 

2. As the concentration of the minority phase is increased the system prefers to phase 
separate into a lamellar and a one-component phase. The lamellae minimise their free 
energy by forming a spiral pattern around the one-component hole. Again this state 
can only be reached by the fluid if hydrodynamic flow is allowed in the system. 

3. When one component of the binary fluid predominates droplets of the minority fluid 
form. Their size is determined by the balance of the surface tension terms in the free 
energy. Phase separation in this system is unaffected by hydrodynamics. 

Sections 2 and 3 of the paper are devoted to a description of the model and the rele- 
vant thermodynamics and to a summary of the numerical approach respectively. Section 4 
summarises phase separation in a 50:50 fluid mixture emphasising the role of hydrodynam- 
ics. The 60:40 mixture which phase separates into a lamellar and one-component state is 
considered in Section 5. Section 6 treats the 90:10 composition ratio where droplets of the 
minority phase are stable. Section 7 summarises the results, putting them in the context of 
previous work, and suggests directions for future research. 



II. MODEL 

We consider a two-dimensional binary fluid with components A, B of number density 
ua, ub respectively. This can be modelled by the Landau free energy! 

* = / dr^f + + + ^(W) 2 } (1) 

where (p = ha — ns is the order parameter of the system and we assume that the total 
density n = + ub is constant. The free energy (JTJ) corresponds to a disordered state at 
high temperatures (a > 0) and an ordered state at low temperatures (a < 0). 



k is related to the surface tension. For k sufficiently large the ordered phase consists 
of two coexisting bulk phases with tia — Ub = ±<£>o> the volume of each determined by the 
initial ratio of ha and As k is decreased the formation of interfaces becomes favourable 
but the fluid remains stable because of the positive curvature energy related to ( in (|1|). For 
f^A ~ n-B the ordered state is then a striped or lamellar phase with the width of the lamellae 
being determined by the interface-interface interaction. For ha <C ub this is replaced by 
a phase of circular droplets of A in B with radius determined by the balance between the 
negative surface tension and the positive curvature free energies of the interfaces. 

Thermodynamic properties of the binary fluid follow directly from the free energy ([[]). 
In particular we shall need the chemical potential Afi which couples to the density difference 

A/2 = ^ = aip + bcp 3 - kV 2 <p + C(V 2 ) V (2) 
dip 

Obtaining the pressure tensor is slightly more complicated^. The pressure parallel to the 
interface follows from ([I]) 

PL = ^ + - MVV) - \ (Vy?) 2 + C^(V 2 ) V - ^(VV) 2 . (3) 

However off-diagonal terms must be added to ensure that the pressure tensor obeys the 
equilibrium condition 

d a P a( s = 0. (4) 

Considering a linear combination of all symmetric tensors having two or four gradient oper- 
ators shows that a suitable choice is 

Pap = {PL + C(VV) 2 + Cd 7 K(VV)}<W 

+Kd a ipd f5 ip - C \d a (pd p {V 2 ip) + dp(pd a (V 2 (p)} . (5) 
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III. LATTICE BOLTZMANN SIMULATIONS 



The aim is to simulate a fluid with equilibrium properties described by the free energy 
(0) which obeys the Navier-Stokes and convection-diffusion equations. To this end we use a 
lattice Boltzmann technique^. 

We consider a square lattice and define two sets of distribution functions {fi} and 
{g{\ on each lattice site x. Each fi,gi is associated with a lattice vector e*j. The re- 
sults presented in this paper are for a 9- velocity model on a square lattice with e^/c = 
(±1, 0), (0, ±1), (±1/V2, ±VV2), (0, 0). 

Physical variables are related to the distribution functions through 

n = fii nU a = (6) 

i i 

^ = Ei9i (7) 

where u is the mean fluid velocity. 

The distribution functions undergo a collision step followed by a streaming step according 
to the evolution equations 

fi{x + e t At,t + At) - fi(x,t) =--(/<- /°), (8) 

gi (x + eiAt,t + At) - gi (x,t) = —{gi-g}) (9) 

where At is the time step, t\ and T2 are relaxation parameters, and ff and are equilibrium 
distribution functions the choice of which determines the physics inherent in the simulation. 
Equations (|8[) and (|9|) are discrete Boltzmann equations with a BGK collision term0. 

Following the standard lattice Boltzmann prescription we assume that ff and can be 
expanded as power series in the bulk velocity 

fi = A + Bu a e ia + Cu 2 + Du a upe ia e i/3 + G a/3 e ia e i( g, (10) 
f° = A + C u 2 , (11) 
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g° = H + Ku a e ia + Ju 2 + Qu a U(3e ia ei(3, (12) 

So = H + J u 2 . (13) 
The expansion coefficients A,B . . . are determined by 



J2f?e ia = nu a , (15) 



/< e ia e i/3 = P a/3 + nu a up, (16) 

Erf = ^ ( 17 ) 



J2di e ia = <PU a , (18) 



E^e^e^ = rA/i<5 Q/3 + <pu a up. (19) 

i 

where P a p and A/x are given by ([|) and (0) respectively and T is a mobility. Note that the 
conditions (0), ([15]), and ( |iT| ) correspond to local conservation of density, momentum and 
density difference respectively. Explicit expressions for the expansion coefficients are given 
in Appendix A. 

Expanding @ and @ to order (At) 2 and using fll4|)-(|l9|) leads to the macroscopic 
equations 

d t n + d a {nu a ) = 0, (20) 
d t {nup) + d a (nu a up) = -d a P a p + i>S7 2 (nup) + c^{A(n)0 a (rai a )}, (21) 
+ a (pu„) = rW 2 A/x - Bd a {^dpPap) , (22) 



where 

9 = (At)c 2 (75 - 1/2) , v = ^-H(At)c 2 , A(n) = (n - ~)At(^ - ^). 

(23) 

It is apparent from these equations of motion that the fluid will evolve through two 
competing growth mechanisms. The first of these is diffusion, described by eqn. ( p2"D with 
u — 0. The second is bulk momentum transfer or hydrodynamic flow described by the 
Navier-Stokes equation (^l|). The relative importance of diffusive and hydrodynamic flow 
can be altered by varying the viscosity v through changes in n (see equation (Pq)). For high 
viscosities the velocities remain sufficiently small that hydrodynamic flow is irrelevant and 
the evolution of the microstructure is diffusive. For low viscosities, however, the Reynolds 
number becomes larger and hydrodynamic effects can dominate in changing the domain 
morphology. Simulations run from the same initial conditions but with high or low viscosities 
enable us to build up a very clear picture of the effects of hydrodynamics on the domain 
growth. 

IV. 50:50 COMPOSITION: METASTABILITY AND THE EFFECT OF 

HYDRODYNAMICS 

We first describe the behaviour of a symmetric binary fluid with a ratio of number 
densities ha '■ ns of 50 : 50 when it is quenched from a disordered state into the ordered 
region. Our aim is to compare the path of the spinodal decomposition for different values of 
the surface tension and of the viscosity. Initial results for this concentration have appeared 
elsewhereS. 

For all runs the system was initialised with n = 1.0 and <p chosen randomly between -0.5 
and 0.5. It was then quenched to a final state defined by parameters a = —1, b = Q = 1. The 
simulations were run with At = 0.004, c = 1, 75 = 0.7786751 and T = 1. The system size 
was 128 x 128 and simulations were typically run for 10 5 time steps. This took approximately 
10 days on a DEC Alpha Workstation. Smaller lattice sizes gave similar results. 
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For k sufficiently negative the equilibrium state is a lamellar phase. After a quench at 
high viscosity the fluid gets stuck in a metastable state of lamellae which have approximately 
the correct width, but which are tangled. For low viscosity, however, hydrodynamic flow 
builds up and can remove topological defects from the system and allow the fluid to reach 
equilibrium. 

Evidence for these conclusions is shown in Figure [l] where snapshots of the domain growth 
are compared for three different sets of parameters: 

(a) high viscosity, two-phase coexistence (n = 50, k = 0.1), 

(b) high viscosity, lamellar equilibrium {t\ = 50, k = —0.85), 

(c) low viscosity, lamellar equilibrium (n = 0.585, k = —0.85). 

Consider first (a) which is at a value of k that corresponds to bulk phase separation 
and a value of T\ that suppresses hydrodynamic flow. After initial transients sharp domains 
form. These are more elongated than for a system with ( = but they grow continuously 
and become more isotropic. An exponent a characterising the growth is commonly defined 
by 

R{t) ~ t a (24) 

where R(t) is a length scale, which we take as the inverse of the first moment of the circularly 
averaged structure factor, and t is the tim At late times a converges to 1/3 as expected 
for Lifshitz-Slyozov diffusive growth in a binary system. Runs at different values of k indicate 
that as k is decreased an increasingly long time is taken to reach the 1/3 regimeS. 

For k < n c ~ —0.8 there is a qualitative change in the behaviour of the system as shown 
in Figure |I]b. These results were obtained for the same value of T\ but with k = —0.85 
where the equilibrium is the lamellar state. Now the system forms portions of lamellae 
of width approximately equal to the equilibrium value. These lamellae are tangled in a 
way that depends on the initial conditions. The diffusive growth to the tangled phase is 
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slow, possibly logarithmic, and once the glassy phase has been reached there is no further 
discernable movement on the time-scale of the simulation. For smaller k the lamellae are 
thinner and the system freezes sooner. 

Remaining with the value k = —0.85 we then ran the same simulation (starting from the 
same initial conditions) for T\ = 0.585. The lower value of r implies a lower viscosity and 
the possibility of hydrodynamic flow at earlier times. The initial behaviour was the same as 
for case (b) with a glassy lamellar phase being formed. However at a later time t there was 
a significant change in the domain pattern. This is most easily seen from the snapshots in 
Figure |l|c where it is apparent that the topological frustration of the glassy phase is being 
removed and the lamellae are lining up in the global equilibrium state. 

An example of the effect of hydrodynamic flow on the defects is shown in Figure 2. The 
lamellae are initially too wide. Hydrodynamic forces tend to lengthen them causing the 
broken lamellar to join and the stripes to buckle. At high viscosities the defect does not 
disappear. 

V. 60:40 COMPOSITION: SEPARATION TO COEXISTING LAMELLAR AND 

BULK PHASES 

We next consider a concentration ratio : % = 60 : 40. Now there is too little of the 
I?-phase to form lamellae of the correct width throughout the system. We shall show that 
the equilibrium state corresponds to phase separation into a lamellar region coexisting with 
a bulk A-phase. As before equilibrium can only be reached with the help of hydrodynamic 
flow. 

Simulations were again run at different values of T\ and k to provide a comparison. 
Snapshots of the time evolution are shown in Figure |3| for the same parameters as those 
used for the symmetric mixture considered in Section IV. Figure 4 is a double logarithmic 
plot of the variation of the domain size with time comparing the growth in the three different 
cases. 



9 



Column (a) of Figure 3 shows the path to bulk phase separation for k = 0.1. Droplets 
form by spinodal decomposition and then grow by Lifshitz-Slyozov diffusion. The curvature 
term in the free energy becomes less important for larger droplets and they become more 
circular. The data in Figure 4 is consistent with the expected growth exponent a = 1/3 at 
late times. 

Figure 3b compares a simulation for k = —0.85, a value for which the lamellar phase is 
stable in the symmetric mixture. The value of t\ is chosen so that only diffusive growth is 
possible. It is apparent from Figure §b that the final state is a mixture of droplets and short 
lamellae. There are no further discernable changes in morphology on the time scale of the 
simulation as confirmed in Figure 4. 

Evidence that this droplet state is metastable is provided when the same simulation is 
run with a low viscosity. Bulk fluid flow now allows the droplets to join to form lamellae 
and align. A surprising amount of movement is seen leading to the growth process shown in 
Figure [3]c. In Figure 4 the onset of hydrodynamic flow in this system is marked by a rather 
sharp increase in the measured length scale. 

The final state is coexistence between a lamellar and a bulk A-phase. The lamellae order 
in a distinctive spiral around the hole of A-phase, to minimise the elastic free energy which 
would be generated by incorrect lamellar spacings. 

Finally we remark that very similar results were obtained for a concentration ratio of 
tia '■ n B = 70 : 30. 

VI. 90:10 COMPOSITION: THE DROPLET PHASE 

Finally we describe phase separation in a highly asymmetric binary fluid with : % = 
90 : 10. Now the equilibrium for k sufficiently negative comprises droplets of B in A. Figure 
5 compares results for two different values of k: 

(a) high viscosity, two-phase coexistence (ri = 50, k = 0.1) 

(b) high viscosity, lamellar equilibrium (n = 50, k = —0.85) 
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The main difference between the growth processes in Figures 5a and 5b is a direct con- 
sequence of the final equilibrium state. For k = 0.1, once the domains have formed they 
continue to grow slowly by the diffusion of material between them. This is the Lifshitz- 
Slyozov growth process, driven by the difference in chemical potential between droplets of 
different size. For k = —0.85 however there is a preferred size for droplets, set by the 
competition between the surface tension and the curvature energy, and growth stops once 
the droplets have reached this size. A second difference is that for the negative value of k 
droplets form much more quickly in the early stages of growth. This is a consequence of the 
reduced surface tension. 

Runs for a low viscosity and negative k showed a negligible influence of hydrodynamics 
on growth for this concentration. This is as expected. Hydrodynamic flow can act to make 
a droplet circular because of the pressure difference between points of different curvature. 
However once the drops are circular hydrodynamic flow cannot directly lead to droplet 
coalescence although it may speed up the diffusive growth^. 



VII. CONCLUSIONS 

We have shown that there is a wealth of interesting behaviour as even rather simple struc- 
tured fluids attain equilibrium. Competition between diffusive and hydrodynamic modes 
lead to final states dependent on the dynamic parameters of the system. In particular we 
have emphasised that hydrodynamic flow is often important in allowing a system to reach 
its equilibrium state. 

The results were obtained using lattice Boltzmann simulations. The approach has two 
particular advantages in this context. Firstly equilibrium is determined by a free energy 
which is an intrinsic part of the simulation so the structure of the equilibrium state can be 
chosen rather naturally. Secondly the viscosity of the fluid can be tuned over a wide range. 
We caution however that the lattice Boltzmann evolution is not described by an H-theorem 
and thus no proof exisits that the path to equilibrium is a physical one. 
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Lattice Boltzmann simulations do not intrinsically include noise although this can be 
added in an ad hoc way. An important question is whether such fluctuations can relax the 
disordered lamellar states. We ran simulations aimed at investigating this but found no 
effect of noise on the glassy structure. In real two-dimensional smectics fluctuations destroy 
the lamellar order. However, for this model we expect a stable lamellar phase, both because 
of the lack of fluctuations and because we impose a mean-field free energy. 

Our conclusions are in broad agreement with those of Bahiana and Oono0 who consid- 
ered the spinodal decomposition of a model of block copolymers designed to give a lamellar 
equilibrium. Despite the long-range interactions in the model a tangled lamellar phase was 
formed after a quench. This could be ordered by including terms approximating hydrody- 
namic flow in the simulation. Other related workEI uses a Langevin approach which includes 
hydrodynamics to model phase separation in microemulsions. The domain growth slows as 
the surfactant density is increased. This method provides an alternative numerical way 
to study complex fluids and it would be interesting to gain a fuller understanding of the 
applicabilities of the different approaches. 

There are many questions that remain to be considered. These include the effect of 
confinement or shear on the phase separation, dynamical asymmetry in the viscosities of 
the two fluid components, and the role of a very low diffusivity, which has been shown to 
alter the path to bulk phase separation in binary fluidsEl. Work is in progress to include a 
surfactant as a third phase rather than modelling its effect by changing the surface tension. 
Extensions to three dimensions are highly desirable. These are feasible but at the limit of 
current numerical resources. 
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APPENDIX A 

A suitable choice of the coefficients in the expansions of the lattice Boltzmann equilibrium 
distributions (10)— (13) consistent with the conditions (14)-(19) is 

A = ! (f V 2 + f V A ~ np (VV) + & (V 2 ) 2 <p + | (VV) 2 ) /(12c 2 ), 
A = n-16A, B = 5n/(12c 2 ), 
C = -2n/(3c 2 ), C = -5n/(24c 2 ), D = 5n/(8c 4 ), 

G XX = -Gyy = f ((d x <p) 2 - (dylf) 2 ) /SC^ + 5( (dylfdy (V^) ~ d x ifd x (V' '(f)) / 8c\ 

G xy = 5nd x (pd y (p - 5( {d x ipd y (V 2 y?) + d y (pd x (V 2 y?)} , 
H = ip-4H, H = 5rA/i/(12c 2 ), K = 5(^/(12c 2 ), 
J = -2v?/(3c 2 ), J = -5vV(24c 2 ), Q = 5^/(8c 2 ). (25) 
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FIGURES 



(a) 



(b) 








FIG. 1. Snapshots of the growth of domains with time for a binary mixture symmetric in 
composition. Grey-scaling from black =>■ white corresponds to 93 = — 1 1. Each column 

represents a different physical situation: (a) a quench to the homogeneous two-phase region k > k c ; 
(b) a quench to the lamellar phase ft < k c in a high viscosity fluid. The lamellae form in a tangled 
pattern which becomes frozen in time; (c) a quench to the same value of k as (b) but for a low 
viscosity fluid. Hydrodynamic modes allow the lamellae to reorder giving, locally, well-defined 
striped regions. 
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FIG. 2. Defect being removed by flow. The lamellae are initially too wide. Hydrodynamic forces 
tend to lengthen them causing the broken lamellar to join and the stripes to buckle. Grey-scaling 
from black =>• white corresponds to 93 = — 1 1. 
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(a) (b) (c) 




FIG. 3. Snapshots of the growth of domains with time for a binary mixture slightly asymmetric 
in composition. Grey-scaling from black =>- white corresponds to ip = — 1 =>- (f = 1. (a) A quench 
to the homogeneous two-phase region k > k c ; (b) a quench to the lamellar phase k < k c in a high 
viscosity fluid. Short lamellae form in a pattern which becomes frozen on the time scale of the 
simulation; (c) a quench to the same value of k as (b) but for a low viscosity fluid. Hydrodynamic 
flow allows the system to attain its equilibrium of coexisting lamellar and bulk- A phases. 
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FIG. 4. Double logarithmic plot of the evolution of the inverse first moment of the structure 
factor as a function of time for each of the simulations shown in Figure 3: ( x ) bulk phase separation, 
(□) a quench to the lamellar phase with high viscosity, (A) a quench to the lamellar phase with 
low viscosity. The straight line has slope 1/3. 
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FIG. 5. Snapshots of the growth of domains with time for a binary mixture highly asymmetric 
in composition. Grey-scaling from black =>- white corresponds to ip = — 1 =>■ </? = 1. (a) A quench 
to the homogeneous two-phase region k > re c ; (b) a quench to the droplet phase k < k c where the 
final length scale is set by the competition between the surface tension and the curvature energy. 
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